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1.1 Introduction 


Many objects in nature are best described geometrically as fractals, with self-similar 
features on all length scales. The universe consists of clusters of galaxies, organized in 
clusters of clusters of galaxies and so on [1.1]. Mountain landscapes have peaks of all 
sizes, from kilometers down to millimeters. River networks consists of streams of all sizes. 
Turbulent fluids have vortices over a wide range of sizes. Earthquakes occur on structures 
of faults ranging from thousands of kilometers to centimeters. Fractals are scale-free in 
the sense that viewing a picture of a part of a fractal one can not deduce its actual size if 


a yardstick is not shown in the same picture. 


If fractals are indeed the geometry of nature, one must still understand how nature 
produces them. A good deal of effort has been put into the geometrical characterization of 
these objects, but there has been practically no progress in understanding their dynamical 
origin. We have a tendency to think of the universe and the crust of the earth as static 
structures because the dynamics that formed these structures were of a much longer time 
scale than the observation period, which could be a human lifetime. The earthquakes that 
we observe last at most a few seconds, whereas the fault formations appear static and were 


built up over millions of years. 


The origin of fractals is a dynamical, not a geometrical, problem. The laws of physics 
are local, but fractals are nevertheless organized over the furthest distances. The mystery 
is enhanced by the fact that large equilibrium systems, operating near their ground state, 
with tend to be only locally correlated. Only at a critical point where a continuous phase 


transition takes place are those systems fractal. 


But real systems are dissipative, that is they have friction, and rarely go to their 
ground state, unlike the ideals discussed in freshman physics. Consider, for example, a 
pendulum. The ideal motion is periodic and for small amplitudes is well approximated 
by a sine wave. To make it more realistic one can put in a drag term, giving rise to 
a damped oscillatory behavior with with a decreasing amplitude theoretically continuing 


forever. However, in the real world the motion will be impeded by some imperfections, 
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perhaps in form of grit. Once the amplitude gets small enough, the pendulum will suddenly 
stop, and this will generally occur at the end of a swing where the velocity is smallest. This 
is not the state of smallest energy, and indeed the probability is a minimum for stopping 
at exactly the bottom of the potential. In a sense, the system is most likely to settle near 


a “minimally-stable” state, far from any “thermal equilibrium.” 


Generalizing to a multi-dimensional system of many coupled pendula, a new issue 
arises. A minimally stable state will be particularly sensitive to small perturbations which 
can “avalanche” through the system. Thus, small disturbances could grow and propagate 
through the system with little resistance despite the damping and impediments. Since 
energy is dissipated through this process, the energy must be replenished for avalanches 
to continue. The systems that we shall study are ones where energy is constantly supplied 


and eventually dissipated in the form of avalanches. 


The canonical metaphoric example is a simple pile of sand. Adding sand slowly to a 
flat pile will result only in some local rearrangement of particles. The individual grains, 
or degrees of freedom, do not interact over large distances. Continuing the process will 
result in the slope increasing to a critical value where an additional grain of sand gives rise 
to avalanches of any size, from a single grain falling up to the full size of the sand pile. 
The pile can no longer be described in terms of many local degrees of freedom, but only a 
holistic description in terms of one sandpile will do. The distribution of avalanches follows 


a power law. 


“Self-Organized Criticality” (SOC) refers to this tendency of large dissipative systems 
to drive themselves to a critical state with a wide range of length and time scales [1.2-4]. 
The idea provides a unifying concept for large scale behavior in systems with many degrees 
of freedom. It has been looked for in such diverse areas as earthquake structure, economics, 


and biological evolution. 


The critical state is an attractor for the dynamics. If a slope were too steep one would 
obtain a large avalanche and a collapse to a flatter and more stable configuration. On 
the other hand, if it were too shallow the new sand will just accumulate to make the pile 
steeper. If the process is modified, for instance by using wet sand in stead of dry sand, the 
pile will modify its slope during a transient period and return to a new critical state. It 
is this resiliency which suggests that self-organized criticality might be a quite universal 


phenomenon. If one builds snow screens locally to prevent avalanches, the pile will again 


respond by locally building up to steeper states, and large avalanches will resume. The 
large fluctuations associated with the avalanches are unavoidable; this might provide some 


food for thought when applied to vast economic or political systems. 


Self-organized criticality complements the concept of “chaos” wherein simple systems 
with a small number of degrees of freedom can display quite complex behavior. Chaos 
is associated with fractal “strange” attractors in the phase space spanned by non-linear 
systems with only a few degrees of freedom. These self-similar structures need have little 
to do with fractals in real spatially extended physical systems. Specifically, chaotic systems 
exhibit white noise with short temporal correlations whereas fractal systems are expected 
to have long range temporal correlations. In contrast, self-organized criticality emphasizes 


unifying features in the coherent evolution of systems with many degrees of freedom. 


1.2 Simulations of Sandpile Models 


The convergence to the self-organized critical state can be demonstrated by computer 
simulations on toy sandpile models. The simplest example is a cellular automaton formu- 
lated on a two-dimensional regular lattice of N sites. Integer variables z; on each site 2 
are used to represent the local sandpile height. Here we consider a two-dimensional lattice 
with open boundaries. Addition of a sand particle to a site 7 is represented by increasing 
the value of z; at that site by unity. When the height somewhere exceeds a critical value 
Zer, here taken to be 3, there is a toppling event wherein 1 grain of sand is transferred 
from the unstable site to each of the 4 neighbor sites, i.e the value of z; is reduced by 4 
and the values of z at the 4 neighbor sites are increased by 1. The updating is done con- 
currently, with all sites updated simultaneously. The initial toppling may initiate a chain 
reaction, where the total amount of topplings is a measure of the size of an avalanche. As 
an example of a state in this model, in Figure 1.0 we show the resulting pile of sand from 
dropping 49, 152 grains of sand on a single site. Already we see signs of a fractal structure 


emerging. 


To explore self organized criticality in this model, one can randomly add sand and 
have the system relax. The result of such an addition becomes unpredictable, with one 


only being able to find the outcome by actually simulating the resulting avalanche. 


In Fig. 1.1 we show a typical state of the sandpile after a large amount of sand has 


been dropped pseudo-randomly. The lattice here is 286 by 184 sites and the boundaries are 
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Fig. 1.0 A sandpile obtained by adding 49,152 grains of sand to a single point. The 


heights z; from 0 through 3 are represented by white, red, blue, and green, respectively. 


open. At first glance, the system appears quite random. There are, however, some subtle 
correlations. For example, never do two adjacent black cells lie adjacent to each other, nor 
does any site have four black neighbors. This follows from the fact that in tumbling a site 


to height 0, a grain of sand is dumped onto each neighbor. 


To a configuration a small amount of sand was added to a site near the center, trig- 
gering an avalanche. We follow this avalanche in Figure 1.2. To trace the avalanch, we 
give the cells which have collapsed a cyan color. Figs. 1.2a and 1.2b show intermediate 
active stages in the collapse. Yellow sites inside or outside the old avalanche area are still 


active. Fig. 1.2c displays the final stable configuration. 


This was a particularly large avalanche, selected for illustrative purposes. The initial 
state was not exactly that of Figure 1.1 because in making this figure we first produced 
several uninterestingly small avalanches. Indeed the ultimate size of the disturbance is 
unpredictable without actually running the simulation. Some cascades may involve a 


single tumbling, and others collapse most of the system. 
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Fig. 1.1 A typical critical state for the sandpile on a 286 by 184 lattice with open 
boundaries. The heights z; from 0 through 3 are represented by black, red, blue, and 


green, respectively. 


Note in Fig. 1.2b the appearance of islands which have not yet collapsed. These have 
all disappeared in Fig. 1.2c, the final relaxed state. This is an exact result, special to this 
model; once in the critical ensemble, it is impossible to trigger a set of avalanches which 


will leave an isolated island of unaltered cells surrounded by disturbed ground. 


Figure 1.3 shows a log-log plot of the distribution of the avalanche sizes s and durations 


t. The linearity indicates a power law, 
T21, (1.1) 


where s is the number of tumblings in an avalanche and P is the probability distribution for 
avalanches of a given size. Thus, the state in Figure 1.1, which at first appears featureless, 
is actually remarkably correlated. For a random distribution of z’s one would expect the 
chain reaction generating the avalanche to be either sub-critical, in which case the avalanche 
would die after a few steps and large avalanches would be exponentially unlikely, or super- 


critical, in which case the avalanche would explode with a collapse of the entire system. 
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Fig. 1.2a-c The progress of an avalanche obtained by increasing one of the z; in Figure 
1.1 to 4. The undisturbed stable sites are colored as in Fig. 1.1 while the still active sites 


in parts a and b are various shades of tan, and the tumbled region is cyan. 


The power law indicates that the reaction is precisely critical, i.e the probability that 
activity at some site branches into more than one active site is balanced by the probability 


that the activity dies. Thus, by evolving through avalanche after avalanche, the matrix 
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Fig. 1.3 A log-log plot of the distribution of the avalanche sizes s and durations ¢ for 
the sandpile model. 


has “learned” to respond critically to the next perturbation. 


Several other quantities which obey fractal scaling laws can be defined for the sandpile. 
For instance, the duration t, that is the number of updatings for an avalanche to complete, 


has a distribution 


P(t)st-™, m2 2.14, (1.2) 
Also, the number of distinct tumbled sites, sg, which is different from the total number of 
topplings since some sites topple more than once, goes as 


P(sa)~ sq "4, Ta & 2.07. (1.3) 


For some applications it is useful to define conditional expectation values. For instance, 
one can define an exponent for the expectation value of the duration, t, for an avalanche 


of size s 


tes", ygs 0:61: (1.4) 


Similar relations connect mutually all the quantities s, t, sg, and the linear size r of 


the clusters [1.5]. The exponent relating sq to T, Ys,r, appears to be 2, indicating that the 
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clusters are compact. The existence of relations between the different stochastic variables 


suggests scaling relations between the exponents. In general 
Ta = 2+ (Ty — 2)/Yay- (1.5) 


As long as the variables x and y are reasonably correlated, this is a mathematical identity 


and does not imply new physics. The compactness of clusters thus suggests (7 — 2) = 
(Tr — 2)/2. 


In order to calculate the power spectrum it turns out that one should calculate the 


expectation value of s? vs duration t, that is 


A(t) = 5° 8? P(s,t). (1.6) 
This quantity is given by the power law 


Mitt, pel. (1.7) 


The values of the exponents quoted here were calculated by Kim Christensen [1.7]. 
The model can be defined in d = 3,4, etc. dimensions. For instance, T ~ 2.31 in three 
dimensions; thus, the values of the exponents depend on d. For an excellent discussion of 


these exponents and their relation to each other, see the paper by Christensen et al. [1.6]. 


1.3. Abelian Sandpile Models 


It would be highly desirable to have an analytical theory, such as the renormaliza- 
tion group theory for equilibrium critical phenomena, by which one could estimate the 
exponents and at the same time gain insights into the mechanisms of self-organized crit- 
icality. We are not yet at that point. However, in a series of papers, Deepak Dhar and 
co-workers have shown that the sand model has some rather remarkable mathematical 
properties [1.8-11]. In particular, the critical attractor of the system is characterized in 
terms of an Abelian group. The properties of the group can be utilized to calculate the 
number of states belonging to the critical attractor, and the rate of convergence to the 
attractor. Further consequences of the Abelian algebra have been explored by one of us 


[1.12-13], and in the following we shall generally follow the discussion given in Ref. [1.13]. 


1.3.1 The Abelian Group 


Dhar introduced the useful toppling matrix A; ; with integer elements representing 
the change in height, z at site į resulting from a toppling at site j [1.8]. Under a toppling 
at site j, the height at site 1 becomes z; — A; j. For the simple two dimensional sand model 


the toppling matrix is given as 


Ai j;=4 i=j 
Aij =—1 i,j nearest neighbors (1.8) 


A; ; =0 otherwise. 


For this discussion there is little special to the specific lattice geometry; indeed the 
following results easily generalize to other lattices and dimensions; in fact, on a Cayley 
tree the model can be solved exactly. The analysis requires only that under a toppling of a 
single site 7, that site has its slope decreased (A; ; > 0), the slope at any other site is either 
increased or unchanged (A; ; < 0,7 Æ i), the total amount of sand in the system does not 


increase È; Ai j > 0), and, finally, that each site can be connected through topplings to 


some location where sand can be lost, such as at a boundary. 


For the specific case in Eq. 1.8, the sum of slopes over all sites is conserved whenever 
a site away from the lattice edge undergoes a toppling. Only at the lattice boundaries can 
sand be lost. Thus the details of this model depend crucially on the boundaries, which we 


take to be open. A toppling at an edge loses one grain of sand and at a corner loses two. 


The actual value of the threshold zr is unimportant to the dynamics. This can be 
changed by simply adding constants to all the z;. Thus without loss of generality we 
consider zr = 3. With this convention, if all z; are initially non-negative they will remain 
so, and we restrict ourselves to states C belonging to that set. The states where all z; 
are positive and less than 4 are called stable; a state that has any z; larger than or equal 
to 4 is called unstable. One conceptually important configuration is the minimally stable 
state C* which has all the heights at the critical value zr. By construction, any addition 


of sand to C* will give an unstable state. 


We now formally define various operators acting on the states C. First, the “adding 
sand” operator a; acting on any C yields the state «;C' where z; = z; + 1 and all other 
z are unchanged. Next, the toppling operator t; transforms C into the state with heights 


2; where zi = zj — A; j- The operator U which updates the lattice one time step is now 
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simply the product of t; over all sites where the slope is unstable, 


uC = |[#c (1.9) 


where p; = 1 if z; > 4; 0 otherwise. Using U repeatedly we can define the relaxation 
operator R. Applied to any state C this corresponds to repeating U until no more 2; 
change. Neither U nor R have any effect on stable states. Finally, we define the avalanche 


operators a; describing the action of adding a grain of sand followed by relaxation 


At this point it is not entirely clear that the operator R exists; that is it might be that 
the updating procedure enters a non-trivial cycle. We now prove that this is impossible. 
First note that a toppling in the interior of the lattice does not change the total amount of 
sand. A toppling on the boundary, however, decreases this sum due to sand falling off the 
edge. Thus, the total sand in the system is a non-increasing quantity. No cycle can have 
toppling at the boundary since this will decrease the sum. Next, the sand on the boundary 
will monotonically increase if there is any toppling one site away. This can not happen in 
a cycle, thus there can be no topplings one site away from the edges. By induction there 
can be no toppling arbitrary distances from the boundary; thus, there can be no cycle, 
and the relaxation operator exists. Note that for a general geometry this result requires 
that every site be eventually connected to an edge where sand can be lost. With periodic 
boundaries no sand would be lost and thus cycles are expected and observed. We call these 
unphysical models “Escher models” after the artist constructing drawings of water flowing 


perpetually downhill and yet circulating in the system. 


It is useful to introduce the concept of recursive states. This set, denoted R, includes 
those stable states which can be reached from any stable state by some addition of sand 
followed by relaxation. As the minimally stable state C* can be obtained from any other 
state by adding just enough sand to each site to make z; equal to three, it belongs to R. 
Thus, one might conveniently define R as the set of states which can be obtained from C* 


by acting with some product of the operators aj. 


It is easily shown that there exist non-recursive, transient states; for instance, no 


recursive state can have two adjacent heights both being zero. One can also show that the 
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self-organized critical ensemble, reached under random addition of sand to the system, has 


equal probability for each state in the recursive set. 


The crucial result of Refs. 1.8-1.11 is that the operators a; acting on stable states 
commute, and that they generate an Abelian group when restricted to recursive states. 
We begin by showing that the operators commute, that is aja;C = aja;C for all C. First 


we express the a’s in terms of toppling and adding operators 


nı n 
ajaj = (fi m) Qi | II m) ajC (1.11) 
k=1 k=nı+1 
where the specific number of topplings nı and n depend on i, 7, and C. Acting on general 
states, the operators t and a all commute because they merely linearly add or subtract 


heights. We therefore can shift a; to the right in this expression: 


ajaj = (i m) [07107163 (1.12) 
k=1 


Now we rearrange the product of topplings. In the non-trivial case that the a-operators 
render either 7 or j (or both) unstable, the product must contain toppling operators cor- 
responding to those unstable sites. We shift those operators to the right. Those operators 


constitute by definition the update operator, U, so we can write 
ajaj = (II tis) UajajC (1.13) 


The factors within the bracket are the remaining t’s. Now, the update operator 
may leave some sites still unstable, and then the product must include further toppling 
operators; working on those sites, we can pull out another factor of the update operator. 
This procedure can be repeated until we have used all the toppling factors and the state 
is stable. Thus, we can identify the operator within the brackets in Eq. (1.12) as the 


relaxation operator R. But aja;C is the same state as ajajC,, so ajaj = aja;C. 


A trivial consequence of this argument is that the total number of tumblings occurring 
in the operations a;a;C and a;a;C are the same. Of course, if a particular site k tumbles 


it can be caused by either addition; the orders of the tumblings may or may not be altered. 


We now prove that the avalanche operators have unique inverses when restricted to 


recursive states; that is, there exists a unique operator a," such that a;(a; ‘C) = C for all 
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CinR. This implies that the operators a; acting on the recursive set generate an Abelian 
group. For any recursive state C we first find another recursive state such that a; acting 


on it gives C, and we then show that this construction is unique. 


We begin by adding a grain of sand to the state C and allowing it to relax. This 
generates a new recursive state a;C’. Now since the state C is by assumption recursive, 
there is some way to add sand to regenerate C from any given state. In particular, there 


is some product P of addition operators a; such that 
C=P al (1.14) 


But the a’s commute, so we have 
C = a; PC (1.15) 


and thus PC is a recursive state on which a; gives C. 


We must still show that this state is unique. To do this consider repeating the above 


process to find a series of states C, satisfying 
(ai)”Cn = C. (1.16) 


Because on our finite system the total number of stable states is finite, this sequence must 
eventually enter a loop. As the original state C can be generated by running around the 
loop in the opposite direction, so C must itself belong to the loop. Calling the length of 


the loop m, we have (a;)"C = C. We now uniquely define a; °C = GG: 


We now have sufficient machinery to count the number of recursive states. As all 
recursive states can be obtained by adding sand to C* , we can write any state C € R in 


the form 
C= (m a) G (1.17) 


Here the integers n; represent the amount of sand to be added at the respective sites. 
However, in general there are several different ways to reach any given state. In particular, 
adding four grains of sand to any one site must force a toppling and is equivalent to adding 


a single grain to each of its neighbors. This can be expressed as the operator statement 


as = || a (1.18) 


jenn 


12 


where the product is over the nearest neighbors to site 3. We can rewrite this equation 
by multiplying by the product of inverse avalanche operators on the nearest neighbors on 


both sides, thus obtaining 
Ay 
I];* == (1.19) 
j 


where E is the identity operator. This allows us to change the powers appearing in Eq. 
(1.17). If we now label states by the vector n = (nı,na,n3.....nn) we see that two such 
states are equivalent if the difference of these vectors is of the form yy 3; Ai; where the 
coefficients 8; are integers. These are the only constraints; if two states can not be related 
by toppling they are independent. Thus any vector n can be translated repeatedly until 
it lies in an N- dimensional hyper-parallelopiped whose base edges are the vectors Aji, 
j = 1,....... N. The vertices of this object have integer coordinates and its volume is the 
number of integer coordinate points inside it. This volume is just the absolute value of the 
determinant of A . Thus the number of recursive states equals the absolute value of the 


determinant of the toppling matrix A. 


For large lattices this determinant can be found easily by Fourier transform. In par- 


ticular, whereas there are 4N stable states, there are only 


exp Nf —— In(4 — 2qz — 2qy) | = (3.2102..)" (1.20) 
(-r,—r) (27) 

recursive states. Thus starting from an arbitrary state and adding sand, the system “self 

organizes” into an exponentially small subset of states forming the attractor of the dynam- 


ics. 


1.3.2 An Isomorphism 


Following Ref. [1.12], we now look into the consequences of stacking sandpiles on top 
of one another. Given stable configurations C and C” with configurations z; and zi, we 
define the state C@C’ to be the state obtained by relaxing the configuration with heights 


zi + zi. Clearly, if either C or C” are recursive states, so is C C”. 


One can show that under the operation @ the recursive states form an Abelian group 
isomorphic to the algebra generated by the a;. First, the addition of a state C with heights 


Zi is equivalent to operating with a product of a; raised to z;, that is 
Bec = (Tai) B. (1.21) 


13 


The operation ® is associative and Abelian because the operators a; are. 


Since any element of a group raised to the order of the group gives the identity, it 


follows that alô! = 1 alt 


i E. This implies the simple formula a; . The analog of this 


for the states is the existence of an inverse state, -C 
-C = (A| - 1) 8&8 C. (1.22) 


Here, n & C means adding n copies of C and relaxing. The state -C has the property that 
for any state BOC 9 (—C) = B. 


The state J = C © (—C) represents the identity and has the property I 6 B = B for 
every recursive state B. The state which is isomorphic to the operator a; is simply a;l. 
The identity state provides a simple way to check if a state, obtained for instance by a 
computer simulation, has reached the attractor, i.e. if a given state is a recursive state: 
A stable state is in R if and only if C ẹ I = C. The proof is simple. By construction, a 


recursive state has this property. On the other hand, since J is recursive, so is C @ I. 


The identity state can be constructed by taking any recursive state, say C* and 
repeatedly add it to itself to use |A| & C = I. However, on any but the smallest lattices, 
[A] is a very large integer. A simpler scheme is given in Ref. [1.12]. Figure 1.4 shows the 
identity state on a 286 by 184 lattice. Note the fractal structure, with features of many 


length scales. 


Here we present another scheme to construct this state. Figure 1.5a-b shows a sequence 
of configurations obtained by pouring sand in from the boundaries: the heights at the edges 
are kept very high, and the inside is initially empty. A variety of fractal structures emerge 
as the interior slowly fills up. Figure 1.5c shows the final stationary state where the 
sand falling in from the boundaries matches that falling off from the updating. Then we 
reverse the boundary conditions and allow the sand to fall back off. Figure 1.5d shows an 
intermediate state of this procedure. When the sand finally stops falling off, we obtain the 


identity state as shown in Figure 1.4. 


1.3.3 A Burning Algorithm and the q=0 Potts Model 


Majumdar and Dhar [1.11] have constructed a simple “burning” algorithm to check 
and enumerate the configurations belonging to the recursive set. The boundary is included 


as a single site, labeled 0. To decide whether a given configuration belongs to R, consider 
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Fig. 1.4 The identity state on a 286 by 184 lattice. The heights are color coded as in 
Fig. 1.1. 


first, at time t = 0, all the sites unburnt, except 0 which is burnt. An unburnt site 2 
is defined to burn at time t =1 if z; exceeds the number of edges connecting 2 to other 
unburnt sites. Those sites burn at t = 2, and so on. The fires start from the edge and burn 
inward. Once some sites have burned, other sites may become burnable since the number 
of edges is reduced. The path of the fire may be indicated by bonds connecting burnt sites. 
If and only if all sites eventually burn under this algorithm does the state belong to the 


recursive set. 


The burning algorithm can also be described in the following way. For a given con- 
figuration, add one particle to each of the edge sites, two particles to the corners, and 
update according to the usual rules. If the original state is recursive, this will generate an 
avalanche under which each site of the system will tumble exactly once. If the state is not 
recursive, some untumbled sites will remain. Figure 1.6 shows such a process underway. 
Here sites which have already burned are shown in cyan, while the remaining sites in the 
center have not yet burned. The small number of sites shown in light tan are the active 


active burning sites. Note the fractal shape of the interphase. 
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Fig. 1.5a-d A procedure for constructing the identity. In parts a and b an empty table 
is having sand poured on from the edges. In part c the entire table is supercritical. In part 
d the boundaries have been opened, and the sand is running back off. The system finally 
relaxes the state shown in Fig. 1.4. In this sequence, heights are color coded as in Fig. 


1.1, with active sites in various shades of tan. 
It is unclear whether the dynamics of the burning algorithm, intended to identify and 
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Fig. 1.6 The burning algorithm being applied to the state in Fig. 1.1. Burnt sites are 


cyan, burning sites are tan, and the remaining sites are colored as in Fig. 1.1. 


count allowed configurations, has anything to do with the dynamics of avalanche formation. 
Majumdar and Dhar studied the scaling of the extension r of the burning cluster vs time 
and found 

tr”, 2 Sy = 5/4. (1.23) 


The spanning tree problem corresponds to a problem in equilibrium statistical me- 
chanics, namely the q = 0 state Potts model. The Potts model is defined in terms of the 


Hamiltonian 


H = q? X 6(0,,0;) (1.24) 
aj 


where the o are q state spins and the summation is over nearest neighbors on a two 
dimensional lattice. Several results about the Potts model are known from an equivalence 
with the two dimensional Coulomb gas model and from conformal field theory. From this 
analogy Majumdar and Dhar [1.10] argue that the height-height correlation (z(r’)z(r’ — 
r)) — (z)? varies as r~?¢ for large separations r, where x; = 2 is the exponent for energy- 


energy correlation function for the Potts model. The height variable in the sandpile model 
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corresponds to the energy density for the Potts model. This again indicates that there are 
strong correlations built into the seemingly random configurations obtained numerically 


(Figure 1.1) for the sandpile. 


1.4. Real Sandpiles and Earthquakes. 
1.4.1 The Dynamics of Sand 


“Sand” represents a “state of matter” to which not much attention has been paid. 
Sand can be shaped into many different forms; that is it can exist in very many stable 
states, almost all of which are out of the flat equilibrium state. Thus, sand contains 
memory: one can write letters in sand. If a heap of sand is perturbed, for instance by 
adding more sand, by tilting the pile, or by shaking it, the system goes from one metastable 
state to another. In some sense this happens by a diffusion process, but this process is 
very different from the process which relaxes a glass of water to equilibrium when shaken. 
The diffusion of sand can stop at any of many stable states, and the process is a threshold 
process, since nothing happens before the perturbation reaches a certain magnitude. It 
might therefore be reasonable to expect real sandpiles to exhibit self-organized critical 


behavior. 


After a couple of false starts with inconclusive results, there are now several exper- 
iments reporting power law distributions of avalanches [1.14-16]. Figure 1.7 shows the 
results of recent experiments by Grumbacher et al. [1.16] They built small heaps on a 
scale, and monitored the distribution of avalanches of particles falling off the edges. The 
figure shows log-log plots of the normalized distribution function of avalanches. The exper- 
iments were performed using iron spheres (triangles) and glass spheres (circles) of the same 
size. In all cases a power law distribution function was found. In a remarkable experiment 
Bretz et al. [1.15] were even able to image the flow of small avalanches not reaching the 


edge and measure their flow and size. 


The threshold dynamics of sand is a paradigm of many processes in nature. Earth- 
quakes occur only when the stress somewhere on the crust of the earth exceeds a critical 
value, and the earthquake takes the crust from one stable state to another. Economic sys- 
tems are driven by threshold processes; the individual agents may change their behavior 
only when conditions reach a certain level. Biological species emerge or die when specific 


conditions in the ecology are fulfilled. Neurons in a neural network fire when the input 
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Fig. 1.7 The avalanche distribution from recent experiments by Grumbacher et al. [1.16] 


from other neurons reaches a threshold level, etc. Indeed, the ideas of SOC have been 
applied to these and many other natural phenomena., including volcanic activity [1.17] 
and solar flares [1.18]. 


1.4.2 Earthquakes and SOC 


The idea of self-organized criticality as applied to earthquakes [1.19- 23], may be 
visualized as follows. Think of the lithosphere of the earth as a collection of tectonic 
plates, being squeezed very, very slowly into each other. At some time in our geological 
history the stresses may have been small, without large ruptures or earthquakes. During 
millions of years, however, the system evolved into a kind of stationary state where the 
build-up of stress is balanced in average by its release during earthquakes. Because of 
the long evolutionary process, the crust has “learned,” by suitably arranging the building 
blocks at hand into a very balanced network of faults, valleys, mountains, oceans and other 
geological structures, to respond critically to any initial rupture. The earthquake can be 


thought of as a critical chain reaction where the process is just barely able to continue. 


For a simple model, consider a two-dimensional lattice of interacting blocks. The 


initial block structure represents a discretization of the space in much the same way as 
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the lattice in lattice gauge theories of particle physics. The block size does not represent 
an intrinsic length scale in the problem. On each block, at site (i,j) acts a force F; j in 
an unspecified general direction of motion in some fault region. In the beginning, F; j 
may assume small random values. The initial state is not important for the long term 
dynamics. Let the force increase uniformly by a infinitesimal small amount per unit time; 
this simulates the slow driving by the tectonic plate motion. Eventually, the force at some 
site (i, j) will exceed a critical threshold value Fo for rupture. The value chosen for Fo 
may be either uniform or random. The initial rupture is simulated by updating the forces 
at the critical site and the sites of the neighbors at (i, j +1) and (i 1, Ĵĵ): 
Fij — 0 


(1.25) 
Fan — Fan + aF; j 


Where nn denotes the nearest neighbors and a is an adjustable parameter. 


These equations represent the transfer of force to the neighbors. This may cause the 
neighbors to be unstable and a chain reaction to take place. This chain reaction is the 
earthquake. The equations are completely deterministic, with no external noise. We are 
not dealing with a noise- driven phenomenon; on the contrary the physics turns out to 
be stable with respect to noise, t.e. noise is irrelevant. When the earthquake stops, the 
system is quiescent until the force at some other location exceeds the critical value and 
a new event is initiated. The process continues again and again. One observes that for 
some time the earthquakes become bigger and bigger. When one is convinced that the 
system has self- organized into a stationary state, one can start measuring the energies of 
subsequent earthquakes as defined by the total number of rupture events following a single 


initial rupture. 


This model was suggested by Olami, Christensen and Feder [1.22] who realized that 
the picture could be directly related to earlier spring-block models. The value of a is 
directly related to the elastic parameters of the crust of the earth. For a = 1/4 the force 
is conserved, i.e. the amount lost on the unstable site equals the total amount gained by 
the 4 neighbor sites. The criticality in this case has been observed to prevail for values of 
a down to 0.05, with only 20% conservation. This came as a surprise since there was then 
a widespread belief that the lack of conservation would spontaneously generate a length 
scale, i. e. a “characteristic earthquake size”. In fact, it seems that criticality occurs 


generically almost independently from the details of the toppling rule. Another model 
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Fig. 1.8 The distribution of earthquakes for a = 0.20 for the model described in the 


text. The straight line yields a power law distribution with an exponent b = 1 — T = 0.8. 


without conservation which appears to exhibit self-organized criticality is the “game of 
life” [1.23], a cellular automaton representing a society of living and dying individuals. 
Perhaps intermittent fluctuations in evolution, like the extinction of the dinosaurs, can be 


seen as manifestations of SOC. 


Figure 1.8 shows the distribution of earthquakes for a = 0.20. The straight line yields 
a power law distribution with an exponent b = 1 — r = 0.8. The straight line indicates 
that the system has self-organized into the critical state. Indeed, real earthquakes exhibit 
a power law distribution, known as the Gutenberg-Richter law. The slope depends on 
the degree of dissipation, (1/4 — a), so there is no universality of the exponent b in the 
non-conservative case. One should not look for unique b-values in nature. Indeed different 


values have been observed in different geographical areas. 


The power law distribution of earthquakes stems from the fractal nature of the SOC 
state, with correlated regions ranging over all length scales; those correlated regions, gen- 
erated by the long term dynamics, are the equivalent of the active faults or fault segments 
in real earthquakes. The fault structure changes on large geological time-scales. More 


realistic long-range SOC models produce faults which topologically look much more like a 


21 


Fig. 1.9 The distribution of A for a = 0.20 for the earthquake model described in the 
text. The slope of the straight line gives u = 0.8. 


real fractal-like arrangement of two-dimensional faults in a three dimensional matrix [1.21]. 


1.5. 1/f Noise 


One-over-f, 1/f, noise can thought of as a signal arising from the superposition of 
avalanches occurring in a self-organized critical state [1.2]. In order to illustrate how this 
works, consider the weighted lifetime distribution of avalanches A(t) as defined in Section 
1.2. Christensen et al. [1.6] showed that if A has a scaling behavior, A(t) ~ t”, then the 


power spectrum S(f) becomes 


S(f) ~ fT, for —1<p<l. (1.26) 


Figure 1.9 shows the distribution of A for œ = 0.20 for the above earthquake model 
[1.24]. The slope of the straight line gives b + 1 = 1.8. 


Figure 1.10 shows a direct measurement of the power spectrum; an exponent ¢ = 1.75 
was found from the slope of the log-log plot, in reasonable agreement with the value 1.8 


expected from the lifetime distribution function defined above. 


In nature, values of the exponent of the 1/f noise in the interval 0.6- 2.0 have been 
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Fig. 1.10 A direct measurement of the power spectrum in the earthquake model. 


reported [1.25]. A value of 1 corresponds to u = 0; this particular value is obtained for 
a = 0.11 in the model. The exponent does depend on the parameters of the model and 


thus is not universal, unlike the conventional exponents for equilibrium critical phenomena. 


1.6 On Forest-Fires and Turbulence 


A liquid driven by imposing a velocity difference v over a length scale L undergoes a 
transition to a turbulent state with vortices exhibiting a large range of length scales. The 
temporal variations of the velocity at a given spot are intermittent, with large and small 
bursts of activity. The energy is dissipated locally within a short length scale known as 


the Kolmogorov length. 


Mandelbrot [1.1] has suggested that in turbulent systems the dissipation of energy 
is confined to a fractal structure with features of all length scales. This behavior can 
be simulated by a simple forest-fire model [1.26]. Distribute randomly a number of trees 
(green dots) and a number of fires (yellow dots). on a two dimensional rectangular lattice. 
Sites can also be empty. Update the system at each time t as follows: (1) grow new trees 
with probability p from sites that are empty at time t — 1; (2) trees that were on fire at 
time t — 1 die (become empty sites) and are removed at time t; (3) a tree that has a fire 


as a nearest neighbor at time t — 1 catches fire at time t. Periodic boundaries are used. 
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Fig. 1.11 A typical state in the dynamics of the forest fire model described in the text. 
Trees are green, fires yellow, and empty sites are black. The periodic lattice is 286 sites by 


178 sites and new trees are born on empty sites with a probability of 1/32 per time step. 


After a while the system evolves to a critical state with fire fronts of all sizes (Figure 1.11). 
Drossel and Schwabl [1.27] have extended the model by adding a small probability f of 
igniting new fires at each time step. In the limit f/p — 0 the ignitions create forest fires 


where the number of trees burned, s, follows a power law, P(s) ~ s!~7, with 7 22. 


The physics of earthquakes can be described in a similar language. The crust in a fault 
region is driven by imposing a force or a strain over a large length L. In the stationary 
state, the energy is dissipated in narrow fault structures forming a fractal set. The spatio- 
temporal correlation functions for the two phenomena are quite similar although the time 
scales are vastly different. In both cases, the energy enters the system uniformly (zero 
wave vector) and leaves the system locally. There are a lot of similarities between the 
earthquake model studied by Olami et al. [1.22] and the forest fire model. The analogy 
has been explored in some detail by Kagan [1.28]. Maybe it is useful to think of self- 


organized criticality and turbulence as one and the same phenomenon. 
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